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Abstract 

The dynamics of a quantum plasma can be described self-consistent ly by the non- 
linear Schrodinger-Poisson system. Here, we consider a multistream model repre- 
senting a statistical mixture of N pure states, each described by a wavef unction. 
The one-stream and two-stream cases are investigated. We derive the dispersion 
relation for the two-stream instability and show that a new, purely quantum, 
branch appears. Numerical simulations of the complete Schrodinger-Poisson sys- 
tem confirm the linear analysis, and provide further results in the strongly non- 
linear regime. The stationary states of the Schrodinger-Poisson system are also 
investigated. These can be viewed as the quantum mechanical counterpart of the 
classical Bernstein-Greene-Kruskal modes, and are described by a set of coupled 
nonlinear differential equations for the electrostatic potential and the stream am- 
plitudes. 
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1 Introduction 



The great degree of miniaturization of today's electronic components is such that 
the de Broglie wavelength of the charge carriers is frequently comparable to the 
dimensions of the system. Hence, quantum mechanical effects (e.g. tunneling) 
are expected to play a central role in the behavior of electronic components to 
be constructed in the next years. In order to describe these quantum effects, it 
is unlikely that classical transport models will be sufficient. Quantum transport 
equations, such as the Schrodinger-Poisson or the Wigner-Poisson systems 
will therefore be a necessary tool in order to understand the basic properties of 
these physical systems. 

In the present paper, we consider a one-dimensional quantum plasma, where the 
electrons are described by a statistical mixture of N pure states, with each wave- 
function ipi obeying the Schrodinger-Poisson system 

dik h 2 d 2 ^ . 

ih— — = — — — echini , i = l, (1) 

dt 2M dx 2 vv ' ' ' K J 

d 2 (h e N 

where 4>(x, t) is the electrostatic potential. Electrons have mass M and charge — e, 
and are globally neutralized by a fixed ion background with density Uq. Finally, 
we assume periodic boundary conditions, with spatial period equal to L. 

The system of Eqs. ©— (§) takes into account diffraction, which is the most promi- 
nent quantum effect, but neglects dissipation, spin and relativistic corrections. 
These effects may be important in more realistic models for small semiconductor 
devices. Nevertheless, it is useful to consider simplified models that capture the 
main features of quantum plasmas. Indeed, Eqs. (0)-(@) are sufficiently rich to 
display a wide variety of behaviors, as it will be seen in the rest of this work. At 
the same time, this model is still amenable to analytic and numerical treatment. 

A physically equivalent approach would consist in considering a Wigner function 
describing the same mixture. The Wigner function approach is a reformulation 
of quantum mechanics in the classical phase space language ||-|ITJ. The price 



to pay for this otherwise appealing formalism, is that Wigner functions can take 
negative values, and cannot therefore be regarded as true probability distributions. 
However, both for the analytical and the numerical treatment of the problems of 
interest in this paper, the Schrodinger-Poisson formalism is more appropriate. This 
is particularly true for the numerical simulations, since the Wigner formalism is 
cast into a two-dimensional phase space, whilst the Schrodinger-Poisson model 
only requires the discretization of the one-dimensional configuration space. Of 
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course, if the number N of streams is large, the numerical cost for the description 
of the system of Eqs. ([!])- (0) is also considerable. Nevertheless, interesting physical 
phenomena (such as instabilities) can take place even with a few streams (N = 1 or 
2), as it will appear in the rest of this work. More subtle effects (Landau damping, 
for instance) would probably require a larger, although hopefully not prohibitive, 
number of streams (N ~ 20 — 30). 

For the analytical study, the hydrodynamic formulation of the Schrodinger-Poisson 
system is particularly convenient, since it makes direct use of macroscopic plasma 
quantities, such as density and average velocity. Moreover, it enables one to per- 
form straightforward perturbation calculations in the same fashion as in the clas- 
sical case. Let us introduce the amplitude A t and the phase Si associated to the 
pure state ipi according to 

ipi = Aiexp(i Si/h) . (3) 

Both Ai and Si are defined as real quantities. The density and the velocity Ui 
of the %— th stream of the plasma are given by 

Introducing Eqs. into Eqs. fll])-© and separating the real and imaginary 

parts of the equations, we find 

On- 

^T + ^~( n i U i) = °> ( 5 ) 

at ox 

dui dui e d(j) h 2 d ( d 2 (y/nl)/dx 2 \ 

~dt +Ui ~dx = Mdx + 2APdx\ ^/h~ /' ^ ' 

d 2 4> e N 
ox 2 e f^i 

Quantum effects are contained in the pressure-like, ^-dependent term in Eq. (|6|). 
If we set h = 0, we simply obtain the classical multistream model introduced by 
Dawson [11|. Therefore, we shall refer to Eqs. @-(0) as the quantum multistream 
model. 

Equations (H)— (0) constitute the mathematical model used in the rest of this work. 
We focus our attention on the one-stream (Sec. 2) and two-stream (Sec. 3) cases, 
the latter being related to the well-known two-stream instability. The relevant 
dispersion relations are derived, and the unstable branches are identified. We 
also investigate the properties of stationary solutions of the Schrodinger-Poisson 
system, which can be viewed as the quantum counterpart of the classical Bernstein- 



Greene-Kruskal (BGK) modes |T2[]. The analytical calculations are checked against 
time-dependent numerical simulations, shown in Sec. 4. Our conclusions are pre- 
sented in Sec. 5. 
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2 One-stream plasma 



In order to fix the basic ideas, we first consider the one-stream case and take N = 1, 
i.e. a single pure quantum state. For brevity, we write n% = n, U\ = u. We obtain 

dn d , s 

« + ^<™» = °- < 8 > 

9it <9w e 90 h 2 d ( d 2 (^fn)jd x 2 \ 

~dt +U ~dx~ = Ml^ + 2M 2 ~~dx~[ 7^ /' 

|^ = -(n-«o). (10) 
A homogeneous zeroth-order solution for Eqs. flSD-fllOD is provided by 

n = n , u = u , = 0, (11) 

where wo is a constant representing the equilibrium velocity of the stream. The 
linear stability of this solution is obtained by Fourier analyzing Eqs. (3)- fll(f) , 

n = n Q + n' exp(i(kx — out)) , (12) 
u = Uq + u' exp(i(kx — uit)) , (13) 
= (f)'exp(i(kx - uot)) (14) 

Retaining only terms up to first order in n', u' and 0' leads to the dispersion 
relation 

(to - ku ) 2 = lu 2 p + h 2 k A /AM 2 , (15) 

where to p = (noe 2 /Me ) 1//2 is the plasma frequency. For uq = 0, the dispersion 
relation given in |0| is recovered, the term kuo merely representing a Doppler shift. 
As the frequency u is always real, there can be neither instability nor damping of 
the wave. 

The classical analog of this system is the "cold plasma" model, which is known 
to sustain nonlinear oscillations when the amplitude of the initial perturbation is 
smaller than a certain value. Beyond that value, the solution becomes singular in 
a finite time, which is a sign that the model is no longer valid. This phenomenon 
corresponds to the breaking of the plasma wave, due to particle overtaking in the 
phase space. Due to the pressure-like term in Eq. ([]), the quantum solution nevers 
becomes singular, as it was shown by computer simulations [[L4J]. 

Let us now turn our attention to the stationary regimes of the system. If all 
quantities only depend on position, then Eqs. ©-© are reduced to 

d [nu) = , (16) 



dx 



du e d<p h 2 d I d 2 (^/n) / dx 2 \ 

U dx = Mdx~ + 21Wdx~{ 7^ / 
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Equations flT6|)-( 



17\) possess the first integrals 



J = nu, (18) 
Mu 2 h 2 ( d 2 {^j)/dx 2 \ 

corresponding to charge and energy conservation. The constant E can be elimi- 
nated by the global shift <p —>■ <fi + E/e, and therefore we assume E = in the rest 
of this section. Eliminating u, introducing A = y/n and using Poisson's equation, 
we obtain 



dx 2 
d 2 4 
dx 2 




(21) 



It can be easily verified that the J = case cannot sustain small-amplitude, 
periodic solutions. Hence, we assume J = uqUq with uq ^ and introduce the 
following rescaling 

x — , J\ — 



u ^/n 

e<f) „ huo. 



MuV H ~ Mu 2 • (22) 

We obtain, in the transformed variables (omitting the stars for simplicity of nota- 
tion), 



d 



2, 



dx 2 



A 2 - 1 , (24) 



a system that only depends on the rescaled parameter H, which is a measure of 
the importance of quantum effects. 



Notice that the classical limit is singular, in the sense that Eq. ( |23|) degenerates 
into an algebraic equation when H = 0. The equation for the electrostatic potential 
then becomes 

dx 2 ~ v/20 1 • ^ bJ 

Equation (|25| ) has a Hamiltonian character, and corresponds to the equation of 
motion of a "particle" moving in a potential V{4>) = <fi — y/2<j). Using this analogy, 
one can see that Eq. ( p5|) has periodic solutions around the equilibrium <fi = 1/2, 
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provided that the initial condition satisfies < <p(x = 0) < 2. The fact that no 
solution exists for sufficiently large values of the potential is easily understood. A 
large potential fluctuation induces a velocity fluctuation, which can drive u(x) far 
from its nominal value uq. If the potential is sufficiently strong, u(x) can even 
vanish, but in that case the relation nu = J = constant implies an infinite density. 
This is the well-known effect of particle overtaking that occurs in the cold plasma 
model. 



Going back to the quantum mechanical case, we have shown that Eqs. (|23|)- 
( plD describe the inhomogeneous QBGK (quantum-BGK) equilibria of the one- 
component quantum plasma. We are faced with the mathematical problem of 
understanding the qualitative properties of the solutions of a coupled, nonlinear 
system of two second order differential equations depending on a parameter. Only 
few analytical results can be obtained, shown in the rest of this section. 



Equations (123]) - ([24]) can be put into Hamiltonian form by using the variables 

A = iA, <j> = <f>/H. (26) 
Notice that the rescaled amplitude A is a purely imaginary quantity. We have 

d 2 A/dx 2 = -d U/dA , d 2 0/dx 2 = -d U/d4> , (27) 
where U = U(A, 0) is the pseudo-potential 

U ^ = ^ l + A2 ^ + ^A2- (28) 

Since the equations of motion are autonomous with respect to the independent 
variable x, the Hamiltonian formulation immediately gives the first integral (sub- 
scripts denote differentiation) 

J = i(^ + £) + &r(4$, (29) 

which is the Hamiltonian function in transformed coordinates. Transforming back 
to the original variables, one obtains the first integral for Eqs. fl2"3"|)-(P4]) 

1 = \(~ A l + H ~ 2 €) + ^(1 - A 2 )^ - ^ . (30) 
According to Liouville- Arnold theorem [ 15 1 , an autonomous two-degree of freedom 



Hamiltonian system is completely integrable if it possesses two first integrals in 
involution and with compact level surfaces. Even if / has not compact level sur- 
faces, a second constant of motion would be a strong indication of integrability of 
the QBGK spatial dynamics. We have tried to find a second constant of motion 
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for Eqs. (p3|)-(|24]) by a variety of methods (geometrical Noether and Lie symme- 
tries, for instance), but without success. However, numerical integrations of Eqs. 
(p3"D-(f24]) for a wide range of values of H, and different initial conditions, strongly 
suggest that, when bounded solutions exist, they are always regular. An additional 
first integral must therefore exist, although its actual expression may be difficult 
to guess. 

It is interesting to perform a linear stability analysis in order to see in what con- 
ditions the system supports small amplitude spatially periodic solutions. Writing 

A = 1 + A' exp(ikx) , 0=1/2 + 0' exp(ikx) , (31) 



and retaining in Eqs. ([23l)-([24]) only terms up to first order in the primed variables, 
we obtain the relation 

H 2 k 4 -4k 2 + 4 = 0. (32) 

Again, we point out the singular character of the classical limit: for H = 0, Eq. (|32}) 
degenerates into a quadratic equation, with solutions k = ±1. The wavenumbers 
being always real, this corresponds to spatially periodic solutions. When H ^ 0, 
we obtain 

^2*2^ (33) 

For H < 1 (semiclassical regime), both wavenumbers are real, and therefore the 
system can sustain spatially periodic oscillations. For H > 1 (strong quantum 
effects), the solutions are spatially unstable, and grow exponentially. For H = 1, 
the spectrum is degenerate, with associated secular terms. The corresponding 
solution is also spatially unstable, growing linearly with x. We conclude that small- 
amplitude stationary solutions of the one-stream Schrodinger-Poisson system can 
only exist in the semiclassical regime, H < 1. 



3 Two-stream plasma 



3.1 Two-stream instability 



We now turn to the more interesting case of two streams, which classically can 
give rise to instability. For this, we consider Eqs. (0)-(0) with N = 2. We first 
linearize around the equilibrium solution n\ = ri2 = no/2, u\ = —112 = uq, = 0, 
where uq 7^ is a nonzero reference velocity, and then Fourier transform both in 
space and time variables. In terms of the dimensionless variables 

O = u/uj p , K = u k/uj p , H = hujp/Mul , (34) 
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the dispersion relation becomes 



Q 4 - f 1 + 2K 2 + — ^— \n 2 -K 2 il — I I 1 - K 2 + — ^— 1=0. (35) 

Consequently, one obtains 

n 2 = i + K 2 + ± 1(1 + 8ff 2 + 4# 2 ff 6 ) 1/2 . (36) 

The solution for Q 2 has two branches, one of which is always positive and gives 
stable oscillations. The other solution is negative (Q 2 < 0) provided that 

(H 2 K 2 -A)(H 2 K 4 -AK 2 + 4) < , (37) 

Notice that instability, when it occurs, always arises through the marginal mode 
(Q = 0). In the classical case, Eq. (|3T| ) yields K 2 < 1, which is the classical 
criterion for the occurrence of the two-stream instability. In the quantum case, 
Eq. flT7|) bifurcates for H = 1. If H > 1, the second factor is always positive, and 
the plasma is unstable if HK < 2. If H < 1, there is instability if either 

< H 2 K 2 < 2 - 2 VI - if 2 , (38) 

or 

2 + 2Vl - if 2 < # 2 ^ 2 < 4 . (39) 

This yields the stability diagram plotted on Fig. 1. The lower instability zone 
is semiclassical, as it represents an extension of the classical instability criterion. 
The upper instability zone, on the other hand, has no classical analog. The two 
zones merge for H — 1. We call Ka, Kb and Kc the wavenumbers at which the 
growth rate vanishes. These are defined by 

H 2 K 2 AB = 2±2Vl -H 2 , (40) 
H 2 K 2 C = 4 . (41) 

It is easy to verify the following property 

K 2 A + K 2 B = K 2 . (42) 

We shall see that these wavelengths are related to the stationary solutions of the 
Schrodinger-Poisson system (quantum-BGK, or QBGK, modes). ^From the pre- 
vious analysis, it appears that quantum mechanics has a destabilizing effect in the 
semiclassical regime, where more modes are unstable compared to the classical 
case. On the other hand, when H > 2, fewer modes turn out to be unstable than 
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in the classical regime. This is the result obtained by Suh, Feix and Bertrand [IE], 
who found that quantum effects are stabilizing for large enough H . 



It would be interesting to know whether the quantum instability is stronger or 
weaker than the classical one. For this, one should search for the maximum (over all 
wavenumbers) growth rate for a fixed value of H, and compare it to the maximum 
classical growth rate. Although the algebra is rather involved, direct inspection of 
the function ~f(K) for different values of H indicates that the maximum classical 
growth rate is always larger than the maximum quantum growth rate (Fig. 2) 
[on this figure, the intersections with the K axis correspond to wavenumbers Ka, 
Kb and Kq as defined in Eqs. (^0|)-(l4l|), for which the growth rate vanishes]. 
Of course, for some wavenumbers, the quantum growth rate may be larger than 
the classical one, as can be seen from Fig. 2. Note that the secondary maximum 
(between wavenumbers Kb and Kq) is considerably smaller than the one between 
K = and Ka , and goes to zero in the classical limit H —>■ 0. 

In the previous discussion, we have ignored the fact that, due to the periodic 
boundary conditions, the momentum variable p = hk is discrete, and can only 
be a multiple of hk , where k = 2ir/L is the fundamental wavenumber. For the 
previous stability calculations, we have assumed an initial equilibrium solution for 
which the density n is spatially uniform and the velocity is equal to ±u . This 
corresponds to wavefunctions of the type 

ip±(x,0) = exp(± iMu x/h) . (43) 

In order to satisfy the periodicity, one must have 

Mu = nhk , n = 1, 2, 3... (44) 

Using the dimensionless variables introduced earlier, this condition becomes 

HK = l/n , (45) 

which implies HKq < 1. This result sets an upper bound on the fundamental 
wavenumber Kq. On the instability diagram of Fig. 1, this means that the upper 
instability region cannot be accessed for Kq, although of course it can for some of 
its harmonics. Computing the intersection of the curve HK = 1 with the lower 
curve on the diagram (K = Ka), we obtain H 2 = 3/4. Therefore, in the region 
H 2 > 3/4, the fundamental wavenumber is always unstable. 

The previous discussion raises the question of the physical meaning of the upper 
instability region in Fig. 1. In particular, we ask whether it is possible to locate 
the fundamental wavenumber in the stable region (between Ka and Kb), and a 
higher harmonic K m = itlKq (where m is an integer) in the unstable upper zone 
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(between Kb and Kc)- The constraints to be satisfied are 



K B < K m <K c , (46) 

K > K A . (47) 

Taking into account Eq. (f§5|), we can write K m = m/(nH), and the previous 
inequalities become 

2 

2 + 2Vl - H 2 < < 4 , (48) 



n 



2 - 2V1 - H 2 < — . (49) 



The square root can be eliminated by summing Eqs. (f48|)-(49). We obtain 

An 2 - 1 < m 2 < An 2 , (50) 

which cannot be satisfied for any pair of integer numbers (n, m) . In summary, it is 
not possible to excite a harmonic in the unstable upper zone of Fig. 1, without also 
exciting the fundamental mode in the lower unstable region. Therefore, at least in 
the semiclassical regime H < 1, we cannot expect to observe a "purely quantum" 
instability. However, as remarked earlier, for H 2 > 3/4 the fundamental mode is 
unstable, as a result of quantum effects. We stress that the above restrictions are 
a result of the periodic boundary conditions, and do not apply to a truly infinite 
plasma, for which momentum space is continuous. Still, periodic conditions can 
be relevant, for instance, to solid state plasmas, where the periodicity is induced 
by the underlying ion lattice. 



3.2 Stationary solutions — QBGK modes 

Let us now consider the stationary states of the two-stream Schrodinger-Poisson 
system (QBGK modes). If all quantities are dependent only on position, Eqs. 
(H)-(|D for N = 2 possess the first integrals 

Ji = n x u x , J 2 = n 2 u 2 , (51) 

Mu\ h 2 d\^)jdx 2 
E t = ~T; e< P- 77TT , (52) 



mux , 


■h 


= n 2 u 2 , 




Mul 


ecj) 


h 2 d 2 (. 


Jni) j dx 2 


2 


2M 




Mu\ 


ecj) 


h 2 d 2 (. 


i/nj) / dx 2 


2 


2M 





E2 = -77^-^-777 V / - • (53) 

2 2M 

We are particularly interested in the case of two symmetric streams, each carrying 
the same current (with opposite sign) and same kinetic energy. Therefore we write 

Ji = -Ji = n u /2 , E 1 = E 2 = Mu 2 /2 (54) 
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where uq 7^ is a nonzero reference velocity. Let ri\ = A\ and n 2 = A%, and 
transform to the dimensionless variables 

x*=uj p x/u , A* 2 = Ax^/^m, (55) 

0* = ( e + Ei)/Mv% , # = hu p /Mul . (56) 



In these variables, the two conservation laws, Eqs. fl52"D-(p3|), and Poisson's equa- 
tion (^) take the form of a six-dimensional dynamical system (we omit the stars 
for simplicity of notation) 

^ - (57) 



d 



2, 



dx 1 



A{ + Ai-l. (59) 



Equations (|57|)-(|59D constitute a coupled, nonlinear system of three second-order 
ordinary differential equations, depending on the control parameter H. They can 
be put into a Hamiltonian form, using a procedure similar to the one employed 
for the one-stream QBGK equations (|2^)-(|24|). This immediately provides a first 
integral, which is the Hamiltonian function itself (its actual expression is rather in- 
volved, and not particularly illuminating). Just as in the one-stream case, we could 
not find any additional constants of the motion, so that we cannot prove rigorously 
that the system is integrable. However, numerical integrations of Eqs. fl57p-(|59"D 
(see end of this section) strongly suggest that the system is indeed integrable, with 
quasi-periodic solutions. 



For the study of small amplitude oscillations, it suffices to expand Eqs. (p7|) — (p9|) 
in the vicinity of the spatially homogeneous equilibrium A\ — 1/ v^2, A 2 = 1/ y/2, 
<fi = 1/2. After Fourier transforming, the following system is obtained, for the 
perturbed quantities A[, <p' 

(4 - H^K^A'i + V20' = 0, i = 1,2 (60) 
V2(A\ + A' 2 )+K 2 <P' = 0, (61) 

where K = kuo/u p is the dimensionless wavenumber. By searching for nontrivial 
solutions, one obtains the relation 

(H 2 K 2 - 4) (H 2 K 4 -AK 2 + A) = 1 (62) 

Notice that this is the same (with an equality sign) as the previously obtained Eq. 
fl57p . Solutions of Eq. ([52]) represent wavenumbers for which both the real and the 
imaginary part of the frequency vanish, and can be considered as the homogeneous 
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limit of generally inhomogeneous stationary states (QBGK modes). If H < 1 there 
are three such solutions, which are the wavenumbers Ka, Kb and Kc defined in 
Eqs. (PUD-flilD . If H > 1, only the solution Kc survives. The other two solutions 
become complex, so that spatially periodic QBGK modes can no longer exist. 

It is also interesting to look for the eigenvectors corresponding to the eigenvalues 
Ka, Kb and K C - From Eqs. (|60|)- (|6lD we can express the perturbed amplitudes 
A\ in terms of the potential. Two cases are possible : (1) if K = Kc = 2/H, 
then one must have 0' = and A[ = —A' 2 ; (2) otherwise, if K = Ka or Kb, we 
have A[ = A' 2 = V2<p'/(H 2 K 2 - 4). The first case is particularly interesting. It 
means that the mode characterized by wavenumber Kc is a "quasi-neutral" mode, 
in the sense that the associated electrostatic potential is zero to first order. Indeed, 
one can see that no plasma parameters (such as the plasma frequency) enter the 
definition of Kc = 2/H. These modes can be actually accessed, for example by 
choosing the fundamental wavenumber Kq = 1/H, which is the largest admissible 
value for K Q [see Eq. (|45|) 1. In this case, the harmonic 2K = 2/H corresponds to 
the quasi-neutral mode. 

Numerical integration of Eqs. (p7[)-(p^D confirms the previous results. For in- 
stance, it was verified that periodic solutions only exist for if < 1. We take 
H = 0.7 and initialize the amplitudes and the potential (at x — 0) with their equi- 
librium value, plus a small perturbation e, i.e. 0(0) = 1/2+60, Ai(0) = (l+6j)/\/2. 
In agreement with the discussion of the previous paragraph, if we choose 6^ = 
and 6i = —62, the wavenumber Kc — 2.857 is linearly excited and thus dominates 
(Fig. 3), while the potential remains very small. On the other hand, if 6i = 62 and 
arbitrary, modes Ka — 1-08 and Kb — 2.645 are linearly excited (Fig. 4). For 
generic perturbations, all three wavenumbers are excited. Of course these results 
are strictly valid only for infinitesimally small perturbations. For moderate values, 
other modes appear (visible on Figs. 3-4, for which e = 0.02), although the linear 
wavenumbers are still dominant. For even larger perturbations, bounded solutions 
no longer exist. 



4 Time-dependent numerical simulations 

A standard numerical technique has been employed in order to integrate the time- 
dependent Schrodinger-Poisson system. Let us write the Schrodinger equation as 

ih^$- =Kip + <$>i/), (63) 
at 

where K is the kinetic part of the Hamiltonian, and $(x,t) is the potential. The 
Hamiltonian is split into these two parts, and each is treated separately. For the 
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potential part, the solution is trivial 

= ?fj n exp{-i<5>At/h) , (64) 

where ip n = ip(nAt), and At is the timestep. For the kinetic part, we use a 
Crank-Nicolson scheme, which is exact to second order in At, 

/ n+l _ in i i 

i* ? At = 2^)" +1 + ■ (65) 

The kinetic operator K is spatially discretized by using the standard centered 
differences formula. The time evolution is obtained by subsequently applying the 
potential and kinetic steps described above. Poisson's equation is solved with a 
Fast Fourier Transform technique just before the potential step (notice that the 
kinetic step does not alter the spatial density, and therefore the potential). The 
resulting numerical scheme is unconditionally stable and second order accurate in 
both space and time variables. Another crucial property of the scheme is that 
it conserves exactly the integral J \ip\ 2 dx (H]]. A typical resolution used in the 



simulation is N x = 512 points, and timestep At = 0.02. 

The initial condition for the simulations is obtained by applying a sudden sinusoidal 
potential to the equilibrium wavefunctions given in Eq. fl4"5p. In dimensionless 
variables, we have 

i ) ±( x i 0) = 2 -1 / 2 exp(± inK x) explieH^ 1 cos(K m x)) , (66) 

where e and K m = itlKq are the amplitude and the wavenumber of the perturba- 
tion, and n,m are integer numbers. We remind that one must have HK = l/n. 

Several simulations have been run in order to compare with the analytical results 
obtained in the linear regime, with excellent agreement between the two. As an 
example, we use the parameters H = 0.25 and n = 5 (K = 0.8), and perturb the 
fundamental mode (m = 1). Figure 5 shows the evolution of two modes of the 
electrostatic potential. The straight line corresponds to the linear growth rate for 
K , as computed from Eq. (0), and closely matches the measured growth rate. 
At saturation, several modes are present. 

The total momentum distribution F(p) is given by the sum of the square modulus 
of the Fourier transform of each wavefunction, with p = hk. As pointed out earlier, 
momentum space is discrete, with Ap = hko (= 0.2 in the above case). For our 
simulations, the total momentum distribution of the unperturbed wavefunctions is 
simply 

HP) = \&(P ~ Mu ) + \5{p + Mu ) , (67) 

where 5 is the Dirac delta function. During the linear phase, the momentum distri- 
bution remains virtually unchanged, whereas at saturation we observe a significant 
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spreading in momentum space, which extends to p ~ ±2.5Muq (Fig. 6). This is 
similar to the behavior of the classical two-stream instability. 

As was shown in the preceding section, one cannot excite a mode in the unstable 
upper region of Fig. 1, without also exciting an unstable mode in the lower zone. 
However, if only the larger wavenumber is initially perturbed, one can hope to 
see it grow with the correct rate before other unstable modes can be excited. In 
order to do so, we take H = 0.9 and n = 5. The fundamental wavenumber is 
thus K = 2/9, and is of course unstable since Ka = 1-18 in this case. We only 
perturb the harmonic 9Ko = 2, which falls within the upper unstable zone, since 
Kb = 1.883 and Ka = 2.222. Of course, several other modes are also unstable 
(namely, those from K to 5K ), but they are not initially perturbed. Figure 7 
shows the evolution of mode 9K , which closely agrees with the result of the linear 
calculation for the growth rate. Several other linearly unstable modes appear at 
a later time, so that at saturation the spectrum is large. This is reflected in the 
momentum distribution, plotted on Fig. 8. 

The results for the fully quantum regime, and particularly for H ~ 1, are more 
surprising. As an example, let us take H = 1, K = 1, and perturb the fundamental 
mode itself. Instability occurs as expected with the correct growth rate (see Fig. 
9). However, instead of saturating at a certain level, the system appears to decay 
with the same rate, and then to grow again. These nonlinear periodic oscillations 
do not damp even for very long times, as it has been checked numerically. Higher 
harmonics, which are linearly stable, are driven by the fundamental mode, and 
show a similar pattern, although at a lower level. The period of the oscillations 
is not universal, and depends on the amplitude of the initial perturbation, which 
confirms that this is indeed a nonlinear effect. A similar, although less pronounced 
behavior, is also observed for H > 1. The momentum distribution is virtually 
unchanged over the entire duration of the simulation. It must be noted that, for 
parameters such that HK = 1, the minimum nonzero value of the momentum is 
IPminl = Muq, and is therefore equal to the momentum of the unperturbed streams. 
In other words, the streams occupy the lowest possible level in momentum space. 
Although we do not have a detailed explanation for this phenomenon, it appears 
to be an example of a completely reversible quantum system, in which the initial 
condition is almost perfectly reconstructed after one period [IS|. This is to be 



compared with the inherently irreversible classical dynamics, for which returning 
to the initial state after saturation is virtually impossible. 

Finally, we have studied the evolution of perturbed QBGK equilibria. It is partic- 
ularly instructive to consider the case where only the Kc = 2/H mode is present. 
This, as detailed in the previous section, is a "quasi-neutral" mode, in the sense 
that the associated electrostatic potential is zero to first order in the perturbation 
parameter. We construct a weakly inhomogeneous equilibrium by using the form 
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of Eq. (^|) for the wavefunctions ipip- Both the amplitude and the phase should 
be sinusoidal functions. For the amplitudes, we have 

A h2 = -^(l±ecos(K c x)) , (68) 

where Kc = mKo is the wavenumber of the QBGK mode; m therefore represents 
the number of density oscillations. The phases Sip are obtained from the definition 
of Eq. (|]), and remembering that, in dimensionless units, UiUi = ±1/2. One 
obtains, to first order in e 

S 1>2 = ±x-2e Kc l sm(K c x) . (69) 

By virtue of the relations nHK = 1 and Kc = mK = 2/H, we obtain that 
m = 2n must be an even number. 

As an example, we take H = 0.5, Kc = 4. With n — 1, we obtain K$ = 2, 
and the number of spatial oscillations is m = 2n = 2. Notice that, in this case, 
the fundamental wavenumber is stable. If we had chosen n = 2, we would have 
had Kq = 1, which is unstable. With this value of H, it is therefore possible to 
construct a stable QBGK mode displaying at most two spatial oscillations. In order 
to have more oscillations, a smaller value of H should be used. For small values of 
the perturbation parameter e, the simulations confirm the linear analytical results, 
and virtually no evolution is observed; besides, the potential fluctuations stay small 
compared to the density ones. For larger values of the perturbation (e = 0.1), 
some temporal variations are observed (Fig. 10), since the wavefunctions given 
by Eqs. d68|)-(|69|) no longer represent an exact stationary state. However, the 
periodic structure is not destroyed, and the potential fluctuations remain an order 
of magnitude smaller than the density fluctuations. These results show that a 
quantum plasma can support almost stationary, quasi-neutral, periodic solutions. 
These display significant density fluctuations for each stream, but small potential 
fluctuations. They have no analog in a classical two-stream plasma. 



5 Conclusion 



In this work, we have introduced a quantum multistream model to describe some 
physical phenomena arising in quantum plasmas. The quantum multistream model 
may be considered as a discrete version of the Wigner-Poisson model, in the same 
sense as the classical multistream model is a discrete form of the Vlasov-Poisson 
system. Indeed, it is well known |19[ that the Wigner-Poisson system is formally 
equivalent to an infinite set of Schrodinger equations, coupled by a scalar potential 
obeying Poisson's equations. However, it is often more appropriate to work with 
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the hydro dynamic formulation of quantum mechanics, since it makes direct use 
of the same physical quantities that are employed in classical physics (density, 
velocity, pressure). Moreover, the stability analysis and perturbation calculations 
become straightforward in the hydrodynamic formulation. On the other hand, 
we have used the Schrodinger representation for the time-dependent simulations, 
since accurate numerical techniques for this equation are well-known from the 
computational literature. 

For the case of the two-stream instability, it has been shown that the dispersion 
relation possesses three branches, one of a semiclassical character and two of a 
purely quantum nature. It is interesting to observe that even the classical branch 
reveals some unexpected features: for small H, quantum effects tend to enhance 
the two-stream instability. More precisely, some classically stable wavenumbers 
are destabilized for sufficiently large values of H < 1. On the other hand, strong 
quantum effect can yield the opposite result: for H > 2, some classically unstable 
wavenumbers become stable. The purely quantum region of the dispersion rela- 
tion (the upper unstable zone of Fig. 1) cannot however be excited without also 
exciting some wavenumber in the lower (semiclassical) region. This means that a 
purely quantum instability cannot be observed for H < 1. Extensive numerical 
simulations, run for different values of the relevant parameters, wholly support the 
analytical results obtained from linear theory. In the fully quantum case, we have 
observed a surprising, yet unexplained, regime of undamped nonlinear oscillations. 
This is a purely quantum effect, which is probably linked to quantum recurrences 



and echoes |T8 



We have also considered the stationary states of the Schrodinger-Poisson system, 
which can be viewed as the quantum analog of the classical BGK modes. Such 
QBGK modes are described by a nonlinear system of coupled second-order differ- 
ential equations, parametrized by the dimensionless Planck's constant H. Such 
system provides an adequate framework for the analysis of QBGK modes, a task 
which would be rather difficult in the Wigner-Poisson formalism. In particular, nu- 
merical simulations have shown that quasi-neutral, spatially periodic, stationary 
states can be created in the two-stream plasma, and can survive over long times. 

Some interesting questions remain to be addressed. In the present work, we have 
considered in detail only the one and two-stream cases. Further investigations 
are needed to explain the properties of the quantum multistream model when the 
number of streams is large. In this case, the Wigner-Poisson system could be a 
more appropriate model, despite its intrinsic mathematical difficulties. However, 
an intermediate number of streams might be sufficient to describe the main physical 
phenomena. Linear and nonlinear Landau damping, for example, should be good 
candidates to test these ideas, both analytically and numerically. 
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Figure Captions 



Figure 1 : Stability diagram for the two-stream plasma. The filled area is unsta- 
ble. The dashed line corresponds to HK = 1. Lower and middle solid curves: K\ 
and K B as defined in Eq. (40). Upper solid curve: K c as defined in Eq. (f4lf). 



Figure 2 : Plot of the squared growth rate 7 as a function of the dimensionless 
wavenumber K, for different values of Planck's constant. H = 0.5, solid line; 
H = 1, dashed line; H = 2, dotted line. The intersections of these curves with the 



K axis correspond to wavenumbers Ka, Kb and Kc as defined in Eqs. (^0|)-(p]). 
For H = 0.5, only the intersection at Ka — 1.035 is shown; the intersections at 
Kb — 3.864 and Kc = 4 are outside the K axis range. For H = 1, the intersections 
are at Ka = Kb = \/2 and Kc = 2. For H = 2, there is only one intersection at 
K c = 1- 

Figure 3 : Stationary solution of the two-stream Schrodinger-Poisson plasma 
(QBGK mode), with H = 0.7, = 0, ei = — e 2 = 0.02. (a) Spatial variation of 
the density fluctuations A[ (solid line), A' 2 (dashed line), and potential fluctuations 
<// (dotted line). Notice that the potential remains small, (b) Fourier transform of 
A[: the linear wavenumber Kc — 2.857 is dominant. 

Figure 4 : Stationary solution of the two-stream Schrodinger-Poisson plasma 
(QBGK mode), with H = 0.7, = 0, t\ = e 2 = 0.02. (a) Spatial variation of the 
density fluctuations A[ (solid line), A' 2 (dashed line), and potential fluctuations 4>' 
(dotted line). Note that the solid and dashed lines are superposed, since A[ ~ A' 2 . 
(b) Fourier transform of A' x : the linear wavenumbers Ka — 1.080 and Kb — 2.645 
are dominant. 

Figure 5 : Two-stream instability - evolution of the fundamental mode Kq = 0.8 
(solid line), and first harmonic 2K (dashed line) of the electrostatic potential, for 
H = 0.25. The straight line corresponds to the growth rate of K computed from 
linear theory, 7 = 0.3116. 

Figure 6 : Momentum distribution for the same case of Fig. 5, at times t = 
(dotted line) and t = 80 (solid line). The final distribution has been magnified by 
a factor two. Momentum space is discrete with Ap = hko = 0.2. 

Figure 7 : Two-stream instability - evolution of mode 9Kq of the electrostatic 
potential, for H = 0.9 and Kq = 2/9. The straight line corresponds to the growth 
rate computed from linear theory, 7 = 0.1085. 

Figure 8 : Momentum distribution for the same case of Fig. 7, at times t = 
(dotted line) and t = 160 (solid line). The final distribution has been magnified 
by a factor two. Momentum space is discrete with Ap = hko — 0.2. 



19 



Figure 9 : Two-stream instability - evolution of the fundamental mode K — 1 
(solid line), and first harmonic 2K (dashed line) of the electrostatic potential, for 
H — 1. The straight line corresponds to the growth rate of mode Kq computed 
from linear theory, 7 = 0.2297. 

Figure 10 : Evolution of a strongly nonlinear (e = 0.1) quasi- neutral (K = Kq = 
4) QBGK equilibrium, for the set of parameters: H = 0.5, K = 2. The number 
of spatial oscillations is m = 2. The density fluctuations n' x (solid line) and n' 2 
(dashed line), and potential fluctuations 0' (magnified by a factor 100, dotted 
line), are shown at different times. 
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